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Abstract 

A method based on the Monte Carlo inversion of the Dirac oper- 
ator on the lattice provides low noise results for the correlations en- 
tering the definition of the heavy meson decay constant in the static 
limit. The method is complementary to the usual method of smeared 
sources, avoids the systematic error arising from optimizing the size 
of the smearing volume and is more efficient for the values of lattice 
parameters that we have explored. 
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The knowledge of leptonic decay constants of heavy mesons is crucial 
for the extraction of the Cabibbo-Kobayashi-Maskawa mixing angles. In par- 
ticular the value of the B meson decay constant constraints the shape of the 
unitarity triangle and sets the size of possible CP violation in B decays. Lat- 
tice calculations are attempting since a few years to provide non perturbative 
estimates of this quantity, but are still confronted with various sources of sys- 
tematic errors. One of them originates from the presence in the dynamics of 
a heavy meson system of two very different mass scales. Reasonable values 
for exploring the B meson dynamics with enough resolution of the heavy 
quark propagation and enough volume for the heavy meson wave function 
are a lattice size of the order of 1.5 fm and a lattice spacing of order of 0.02 
fm which lead to about 75 lattice points, a value beyond present computer 
capabilities. 

One reaches the physical region of the heavy quark mass by interpolat- 
ing between the results which are obtained in the charm quark mass region 
and those deriving from an expansion of the fermion action in the inverse of 
the heavy quark mass. The first term in the latter expansion is the so called 
static limit [1] and corresponds to the approximation where the heavy quark 
does not propagate in space. An accurate knowledge of the decay constant in 
this limit is an essential ingredient of the calculation of the physical value. In 
spite of the great computer effort dedicated to its estimate the values quoted 
in recent papers are still affected by a 10 — 15% error, an important fraction 
of which is due to the large statistical fluctuations. In this letter we propose 
an alternative way of measuring the correlator of the heavy meson in the 
static limit which makes possible to reduce the statistical errors significantly. 

On the lattice and in the static limit, the correlator of two local heavy- 
light pseudoscalar bilinears is given by: 

G(t) = j^Y,tyty^np(*,y)^si( X ,y)]) (i) 

where S q is the light quark propagator, P(x,y) is the Polyakov line from the 
point x to y, and L and T are the space and time sizes respectively. 

At large times there is a single state dominating the spectral decompo- 
sition of the correlation function and one can extract the matrix element of 
the heavy-light current between such a state and the vacuum: 



2 



G {t) tl J^ e Z 2 L e- AEt (2) 
The decay constant in the static limit can be derived from Z^: 

jpsjMp = V2Z Ren Z L a~ 3 / 2 (3) 

where Z Ren is the suitable renormalization constant of the lattice current. 

With standard methods the signal of the local currents is too noisy to 
wait for large times where the lowest lying state dominates : the extraction 
of the matrix element is then performed in two steps. In the first one uses 
smeared currents which optimize the projection on the lowest mass state 
and allow an early estimate of the energy shift in eq. 2, in the second, one 
analyses the ratio of single and double smeared correlation functions to get 
the estimate of the local matrix element. Error propagation is minimized by 
performing fits of these correlations globally. A typical source of systematic 
error arises from the size of the smearing volume and from the choice of the 
window in time used for the fits. 

In the static approximation, the final space point of the correlator coin- 
cides with the initial space point. The inversion of the Dirac operator is very 
expensive in computer time and is normally performed for a single initial 
point on each gauge configuration. The major source of statistical errors can 
be ascribed to the difficulty of estimating the correlator in eq. 1 for many 
independent initial points. 

The method presented in this letter is based on a Monte Carlo inver- 
sion of the Dirac operator which allows to perform the average of the static 
correlator over all the initial points, to obtain a visible signal at large enough 
times and to avoid the use of smeared sources. We introduce an auxiliary 
action with pseudofermion fields: 

W^] = £IKM(*)I 2 ( 4 ) 

X 

where (j)(x) is the pseudofermion field, and for the lattice Dirac operator D 
we follow the standard Wilson formulation: 

1 1 3 

[Q<f)](x) = nf 6 D</> = — 75<f>(x) ~ -75 £ U »( X )( l ~ 1/M X + A*) 

~7bE^-A*)(1+7m)^-m) (5) 
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where k, is the Wilson hopping parameter related to the bare mass. 

The sum ^[^7] + S$\U,4>\, where S g is the standard Wilson action for 
the gauge sector, is the "bermion" action used in the study of the dynamical 
flavour dependence of lattice QCD by extrapolating from negative flavour 
numbers [2]. In the context of this letter, the bermion fields are not dynam- 
ical, they are thermalized in a fixed gauge confuguration and only used to 
estimate fermion propagators, like in the fermion algorithm of ref. [3]. The 
presence of the square of the Dirac operator is needed for the Montecarlo 
update of the pseudofermion fields and implies that the two point function 
of these fields represents the inverse of such a square. The usual inverse 
Dirac operator can be obtained by the remultiplying the pseudofermion field 
correlations with the Dirac operator. We define two correlators, the first 
corresponding to eq. 1: 



G(t) = ^E^tM[Q<f>]\^sP(x,y)^<f>(y)) (6) 

and a second one corresponding to propagators with the square of the Dirac 
operator: 

G {Q2 \t) = ^E4AM t (^(^)^,2/)^)) (7) 

The statistical fluctuations of these operators can be reduced by a "one 
<f> integral" technique analogous to the "one link integral" method [5]: one 
replaces the pseudofermion field by its average in the surrounding pseud- 
ofermion and gauge configuration. The explicit expression for such a re- 
placement is: 

^{x)^ -A-^Q 2 ^- Act>]{x) (8) 

where 



* = ^F- (») 



1 + 16k 2 
Tit 2 

The action for the fields which are replaced in the same observable 
through eq. 8 should be separable, which implies that the improved corre- 
lation functions G(t) and G^ Q '\t) can only be calculated at time distances 
greater than two. 
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As already noticed earlier [4] the correlator constructed from the inverse 
of the square of the Dirac operator in eq. 7 projects more precociously on the 
lowest lying state. We have therefore extracted the energy shift AE from this 
correlator and then fixed its value in the fit of the canonical correlator G{t). 
We have tested this method at two values of (3: 5.7 and 6.0 on a 16 3 X 32 lat- 
tice by measuring correlation functions in all four directions. The simulations 
were performed on a 25 Gigaflop machine of the APE series. The update pro- 
cedure was for the gauge sector a Cabibbo-Marinari pseudo-heatbath [6] fol- 
lowed by three overrelaxation sweeps and for the pseudofermions a heat bath 
followed by ten overrelaxation sweeps [7] . The gauge field configurations were 
separated by 1000 sweeps and for each fixed gauge configuration the calcula- 
tion of the pseudofermion observables was performed every 5 pseudofermion 
sweeps in a total of 2000 pseudofermion updates, i. e. 400 times. In addition, 
for each gauge configuration, pseudofermion measurements were taken after 
500 thermalization updates which were tested to be sufficient to equilibrate 
the pseudofermion system. The Montecarlo inversion of the Dirac operator 
with 400 samplings achieves a modest precision with respect to the inversion 
by a minimization algorithm, but becomes of comparable quality when one 
uses the sum over all space time points. For example, the zero momentum 
pion correlation, averaged over 30 gauge configurations, at k, = 0.165 has 
a 10% error at t = 13 if calculated with the Dirac deterministic inversion, 
and the same error at t = 10 if calculated with the Monte Carlo inversion. 
The sum over all space time points in the case of fs provides the additional 
advantage of sampling also the heavy quark propagator through the whole 
lattice. The computer time required with the 2000 updates and 500 thermal- 
ization sweeps of the pseudofermion system is comparable to that of a value 
of a single origin Dirac inversion at rn^/m 2 p ~ 0.4. 

The results of our simulation refer to a set of 30 gauge configurations. 
The errors are evaluated as follows: we form clusters of the 400 measurements 
of the correlation functions at fixed gauge configurations or larger clusters 
containing 3 gauge configurations, i. e. of a total of 400 X 3 measurements. 
This subdivides the sample of 30 X 400 pseudofermion measurements into 30 
or 10 clusters respectively. On these clusters we apply a jacknife algorithm 
to estimate the errors of the correlations. Each "jacknife cluster", defined by 
single elimination of a standard cluster from the total sample is fitted with a 
MINUIT minimization program. From the spread of the fits on each jacknife 
cluster we extract the error that we quote for our results. The procedure 
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with standard clusters gives consistent results within the quoted errors. 

The results for the effective mass in the case (3 = 5.7 at a value of 
the hopping parameter k, = 0.165 corresponding to a ratio ra^/ra^ of 0.5 are 
presented in figure 1 and compared with the fit. For the operator of eq. 7 a 
one mass fit is generally adequate. For the operator of eq. 6 a two mass fit is 
necessary with the lowest mass value fixed by the fit of the other operator. 
The results for different values of k, and for different choices of the starting 
value of the time window used in the fit are given in Table 1. In the case 
of the single mass fit the window is shifted to larger times to reach a better 
projection on the lowest mass state. 

For Zl we give separately the statistical error of the G(t) correlation fit 
and the error deriving from the uncertainty by which the correlation function 
G& 2 \t) fixes the lowest energy shift value. Our results can be compared 
with those of [8], obtained on a smaller volume (12 3 X 24). Our values, less 
affected by finite size effects, lie in general below, exhibit smaller errors and 
are obtained with a third of the number of gauge configurations used in [8]. 

For the results at (3 = 6.0 one needs a two mass fit also for the corre- 
lation function G^ Q2 \t) as it can be seen from figure 2 which presents the 
effective mass at a value of the hopping parameter k corresponding to a ratio 
m\lm 2 p of 0.5. The continuous line is our two mass fit. Indeed, at (3 = 6.0 
the lattice spacing is smaller and at fixed time distance in lattice units the 
correlation functions are sampled at smaller physical units where the excited 
states are still important. The compensating effect of a smaller energy shift 
in lattice units which would allow the signal to survive up to larger lattice 
times is partly lost because of the linear divergences which contribute with a 
constant to the mass in lattice units. This effect is stronger in the correlation 
of the operator of eq. 6: in this case, the results for Zl, of the two mass fit 
show a systematic decrease with increasing values of the starting time of the 
window used for the fit. This does not happen for the results at (3 = 5.7 
which are insensitive to the choice of the fitting window. The final value 
we quote is obtained from the last time window, but should still be consid- 
ered as an upper bound of the true value. The new preliminary results of 
the APE [9] and MILC [10] collaborations at (3 = 6.0 are indeed compatible 
with ours but lie below within one standard deviation. The extrapolation of 
our results to the chiral limit interpolates rather well between the results at 
(3 = 6.1 and (3 = 5.9 of [8]. A work applying this method to the study of the 
dynamical flavours dependence of fs where the pseudofermions are replaced 
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by the dynamical bermion fields is in progress [11]. 

Our method appears to reach a better statistical precision with respect 
to standard method. The gain in statistics comes from three sources: the 
Montecarlo Dirac inversion which provides an alternative way of measuring 
correlations where, as in the case of static pseudoscalar decay constant, the 
additional average over all points is important, the measurement of a corre- 
lator which leads to an earlier determination of the lowest lying state and 
the "one </>" integral which reduces the statistical fluctuations. The method 
discussed in this paper is complementary to the existing smearing procedures 
and can be conjugated with them to further reduce the statistical and sys- 
tematic uncertainties of lattice determination of static pseudoscalar decay 
constants. 

Acknowledgments. We thank C. Allton and C. Bernard for commu- 
nicating to us their preliminary results prior to publication and A. Vladikas 
for stimulating discussions. 
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A£(* min = 4) 


A£(* min = 5) 


A£(* min = 6) 


AE 


0.163 
0.165 
0.16625 


0.8145(6) 
0.7848(7) 
0.7645(9) 


0.8126(7) 
0.7853(9) 
0.7671(12) 


0.8118(10) 
0.7865(12) 
0.7701(20) 


0.812(2) 
0.786(2) 
0.770(4) 




■^L(^min = 3) 


Z L (t m in = 4) 


■^L(^min = 5) 


Z L 


0.163 
0.165 
0.16625 


0.649(4) 
0.590(7) 
0.558(7) 


0.655(8) 
0.591(10) 
0.559(9) 


0.655(20) 
0.590(20) 
0.560(20) 


0.655(8)[7] 
0.590(7)[7] 
0.559(9)[12] 



Table 1: Results for AE and Zl calculated on a 16 3 x 32 lattice at j3 = 5.7. 
In the last column we quote our final result: for Zl we separate the statistical 
error () from the one coming from the determination of AE []. 





A£(* min = 3) 


A£(* min = 4) 


AE(^ min = 5) 


AE 


0.154 
0.155 
0.1557 


0.635(7) 
0.616(7) 
0.598(10) 


0.635(10) 
0.616(10) 
0.598(14) 


0.632(30) 
0.615(18) 
0.598(20) 


0.635(7) 
0.616(7) 
0.598(10) 




■^L^min = 3) 


Z L (t min = 4) 


■ZL(£ m in = 5) 


Z L 


0.154 
0.155 
0.1557 


0.250(4) 
0.231(3) 
0.213(3) 


0.244(5) 
0.223(5) 
0.207(5) 


0.237(8) 
0.220(7) 
0.200(7) 


0.237(8)[10] 
0.220(7)[7] 
0.200(7)[10] 



Table 2: The same as in table 1 for (3 = 6.0. 
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Figure 1: The effective mass for the correlation function G^ 2 \t) (diamonds) 
is compared with the single mass fit and the effective mass for the correlation 
function G(t) (crosses) is compared with the two mass fit for (3 = 5.7 and 
k = 0.165. 
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Figure 2: The effective mass for the correlation G^ 2 \t) for (3 
k, = 0.155 is compared with the two mass fit. 
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